/**
 * @file   Matrix.h
 * @author xingyifan <xingyifan@xingyifan-VirtualBox>
 * @date   Tue Jul 13 10:43:16 2021
 * 
 * @brief  some Algebra function
 * 
 * 
 */

#include "MyExcept.h"

#ifndef __MATH__
#define __MATH__

#include <cmath>
#include <ctime>
#include <cstdlib>

#endif

using namespace std;

#ifndef __ALGEBRA__MATRIX__
#define __ALGEBRA__MATRIX__

namespace Algebra {
    /**
     * @file   Matrix.h
     * @author xingyifan <xingyifan@xingyifan-VirtualBox>
     * @date   Wed Jul 14 10:34:03 2021
     * 
     * @brief  vector
     * 
     * 
     */

    template <typename T>
    class Vector {
    private:
	int len;
        T * value;
    public:
	Vector() {
	    len = 0;
	    value = NULL;
	}
	Vector(int length) {
	    len = length;
	    value = new T[len];
	    for (int i=0;i<len;i++) {
		value[i] = 0;
	    }
	}
	Vector(int length, T a[]) {
	    len = length;
	    value = new T[len];
	    for (int i=0;i<len;i++) {
		value[i] = a[i];
	    }
	}
	Vector(const Vector& v) {
	    len = v.get_len();
	    value = new T[len];
	    for (int i=0;i<len;i++) {
		value[i] = v[i];
	    }
	}
	
        T &operator[] (int i) {
	    try {
		if (i>=0 && i<len) {
		    return value[i];
		} else {
		    throw MyExcept(ARRAY_BOUND);
		}		
	    }
	    catch (MyExcept &e) {
		e.what();
		return value[0];
	    }
	}
        T &operator[] (int i) const {
	    try {
		if (i>=0 && i<len) {
		    return value[i];
		} else {
		    throw MyExcept(ARRAY_BOUND);
		}		
	    }
	    catch (const MyExcept &e) {
		e.what();
		return value[0];
	    }
	}
	
        Vector operator+ (T delta) {
	    Vector v(len);
	    for (int i=0;i<len;i++) {
		v[i] = value[i] + delta;
	    }
	    return v;
	}
        Vector operator+ (const Vector& _v) {
	    Vector v(len);
	    for (int i=0;i<len;i++) {
		v[i] = value[i] + _v[i];
	    }
	    return v;
	}
	
        Vector operator- (T delta) {
	    Vector v(len);
	    for (int i=0;i<len;i++) {
		v[i] = value[i] - delta;
	    }
	    return v;
	}
        Vector operator- (const Vector& _v) {
	    Vector v(len);
	    for (int i=0;i<len;i++) {
		v[i] = value[i] - _v[i];
	    }
	    return v;
	}

        Vector operator* (T delta) {
	    Vector v(len);
	    for (int i=0;i<len;i++) {
		v[i] = value[i] * delta;
	    }
	    return v;
	 }
	 double operator* (const Vector& v) {
	     T sum = 0;
	     for (int i=0;i<len;i++) {
		 sum += value[i] * v[i];
	     }
	     return sum;
	 }

	 Vector operator/ (double delta) {
	    Vector v(len);
	    for (int i=0;i<len;i++) {
		v[i] = value[i] / delta;
	    }
	    return v;
	 }

	 void operator= (const Vector& v) {
	     if (value!=NULL) {
		 delete [] value;
	     }	     
	     len = v.get_len();
	     value = new T[len];
	     for (int i=0;i<len;i++) {
		 value[i] = v[i];
	     }
	 }
	 bool operator== (const Vector& v) {
	     if (len!=v.get_len()) {
		 return false;
	     }
	     for (int i=0;i<len;i++) {
		 if (value[i]!=v[i]) {
		     return false;
		 }
	     }
	     return true;
	 }
	
	friend ostream& operator<<(ostream& out, Vector& v) {
	    out << "[";
	    for (int i=0;i<v.get_len();i++) {
		out << v[i] << " ";
	    }
	    out << "]";
	    return out;
	}
	
	int get_len() const {
	    return len;
	}
	void set_len(int length) {
	    try {
		if (len==0) {
		    len = length;
		    value = new T[len];
		    for (int i=0;i<len;i++) {
			value[i] = 0;
		    }
		} else {
		    throw MyExcept(VECTOR_LENGTH);
		}
	    }
	    catch (const MyExcept& e) {
		e.what();
	    }
	}

	double one_norm() {
	    double sum = 0;
	    for (int i=0;i<len;i++) {
		sum += fabs(value[i]);
	    }
	    return sum;
	}

	double inf_norm() {
	    double max = 0;
	    for (int i=0;i<len;i++) {
		double tmp = fabs(value[i]);
		if (tmp>max) {
		    max = tmp;
		}
	    }
	    return max;
	}

	double two_norm() {
	    double sum = 0;
	    for (int i=0;i<len;i++) {
		sum += value[i] * value[i];
	    }
	    return sqrt(sum);
	}
	
	~Vector() {
	    delete [] value;
	}
    };

    /**
     * @file   Matrix.h
     * @author xingyifan <xingyifan@xingyifan-VirtualBox>
     * @date   Wed Jul 14 10:34:19 2021
     * 
     * @brief  normal matrix
     * 
     * 
     */

    template <typename T>
    class Matrix {
    private:
	int row;
	int column;
	Vector<T> * value;
	
    public:
	Matrix(int _row, int _column) {
	    row = _row;
	    column = _column;
	    value = new Vector<T>[row];
	    for (int i=0;i<row;i++) {
		value[i].set_len(column);
		for (int j=0;j<column;j++) {
		    value[i][j] = 0;
		}
	    }
	}
	Matrix(int _row, int _column, double A[]) {
	    row = _row;
	    column = _column;
	    value = new Vector<T>[row];
	    for (int i=0;i<row;i++) {
		value[i].set_len(column);
		for (int j=0;j<column;j++) {
		    value[i][j] = A[i*column+j];
		}
	    }
	}
        Matrix(const Matrix& m) {
	    row = m.get_row();
	    column = m.get_column();
	    value = new Vector<T>[column];
	    for (int i=0;i<row;i++) {
		value[i] = m[i];
	    }
	}
	
        Vector<T> &operator[] (int i) {
	    try {
		if (i>=0 && i<row) {
		    return value[i];
		} else {
		    throw MyExcept(ARRAY_BOUND);
		}		
	    }
	    catch (MyExcept &e) {
		e.what();
		return value[0];
	    }
	}
        Vector<T> &operator[] (int i) const {
	    try {
		if (i>=0 && i<row) {
		    return value[i];
		} else {
		    throw MyExcept(ARRAY_BOUND);
		}		
	    }
	    catch (const MyExcept &e) {
		e.what();
		return value[0];
	    }
	}
	
        Matrix<T> operator+ (T delta) {
	    Matrix<T> m(row, column);
	    for (int i=0;i<row;i++) {
		m[i] = value[i] + delta;
	    }
	    return m;
	}
        Matrix<T> operator+ (const Matrix<T>& _m) {
	    Matrix<T> m(row, column);
	    for (int i=0;i<row;i++) {
		m[i] = value[i] + m[i];
	    }
	    return m;
	}
	
        Matrix<T> operator- (T delta) {
	    Matrix<T> m(row, column);
	    for (int i=0;i<row;i++) {
		m[i] = value[i] - delta;
	    }
	    return m;
	}
        Matrix<T> operator- (const Matrix<T>& _m) {
	    Matrix<T> m(row, column);
	    for (int i=0;i<row;i++) {
		m[i] = value[i] - m[i];
	    }
	    return m;
	}

	 Matrix<T> operator* (T delta) {
	     Matrix<T> m(row, column);
	    for (int i=0;i<row;i++) {
		m[i] = value[i] * delta;
	    }
	    return m;
	 }
	 Vector<T> operator* (const Vector<T>& _v) {
	     Vector<T> v(row);
	     T sum;
	     for (int i=0;i<row;i++) {
		 sum = 0;
		 for (int j=0;j<column;j++) {
		     sum += value[i][j] * _v[j];
		 }
		 v[i] = sum;
	     }
	     return v;
	 }
	 Matrix<T> operator* (const Matrix<T>& _m) {
	     Matrix<T> m(row, _m.get_column());
	     T sum;
	     int _row = _m.get_row();
	     int _column = _m.get_column();
	     for (int i=0;i<row;i++) {
		 for (int j=0;j<_column;j++) {
		     sum = 0;
		     for (int k=0;k<_row;k++) {
			 sum += value[i][k] * _m[k][j];
		     }
		     m[i][j] = sum;
		 }
	     }
	     return m;
	 }

	 Matrix<T> operator/ (T delta) {
	     Matrix<T> m(row, column);
	    for (int i=0;i<row;i++) {
		m[i] = value[i] / delta;
	    }
	    return m;
	 }
	 
	 void operator= (const Matrix<T>& m) {
	     if (value!=NULL) {
		 delete [] value;
	     }	     
	     row = m.row;
	     value = new Vector<T>[row];
	     for (int i=0;i<row;i++) {
		 value[i] = m[i];
	     }
	 }
	 bool operator== (const Matrix<T>& m) {
	     if (row!=m.get_row() || column!=m.get_column()) {
		 return false;
	     }
	     for (int i=0;i<row;i++) {
		 if (!(value[i]==m[i])) {
		     return false;
		 }
	     }
	     return true;
	 }
	
	friend ostream& operator<<(ostream& out, Matrix<T>& m) {
	    out << "[" << endl;
	    for (int i=0;i<m.get_row();i++) {
		out << m[i] << endl;
	    }
	    out << "]" << endl;
	    return out;
	}       
	
	double one_norm() {
	    double sum, max = 0;	    
	    for (int i=0;i<column;i++) {
		sum = 0;
		for (int j=0;j<row;j++) {
		    sum += fabs(value[j][i]);
		}
		if (sum>max) {
		    max = sum;
		}
	    }
	    return max;
	}

	double inf_norm() {
	    double max = 0;
	    for (int i=0;i<column;i++) {
		double tmp = value[i].one_norm();
		if (tmp>max) {
		    max = tmp;
		}
	    }
	    return max;
	}

	Matrix<T> transpose() {
	    Matrix<T> m(column, row);
	    for (int i=0;i<column;i++) {
		for (int j=0;j<row;j++) {
		    m[i][j] = value[j][i];
		}
	    }
	    return m;
	}

	double normE() {
	    try {
		if (row==column) {
		    double sum = 0;
		    for (int i=0;i<row;i++) {
			for (int j=0;j<row;j++) {
			    if (i!=j) {
				sum += value[i][j] * value[i][j];
			    }			    
			}
		    }
		    return sqrt(sum);
		} else {
		    throw MyExcept(NOT_A_SQUARE);
		}
	    }
	    catch (MyExcept& e) {
		e.what();
		return -1;
	    }
	}

        double normF() {
	    try {
		if (row==column) {
		    double sum = 0;
		    for (int i=0;i<row;i++) {
			for (int j=0;j<row;j++) {
			    sum += value[i][j] * value[i][j];	    
			}
		    }
		    return sqrt(sum);
		} else {
		    throw MyExcept(NOT_A_SQUARE);
		}
	    }
	    catch (MyExcept& e) {
		e.what();
		return -1;
	    }
	}

	int get_row() const {
	    return row;
	}
	int get_column() const {
	    return column;
	}
	
	~Matrix() {
	    delete [] value;
	}
    };

    /**
     * @file   Matrix.h
     * @author xingyifan <xingyifan@xingyifan-VirtualBox>
     * @date   Wed Jul 14 10:33:43 2021
     * 
     * @brief  sparse matrix
     * 
     * 
     */

    template <typename T>
    class SpMatrix {
    private:
	int l_i, l_j, nnz;       
	Vector<T> * IA;
	Vector<T> * JA;
	Vector<T> * SA;

    public:
	SpMatrix(int _l_i, int _l_j, int _nnz) {
	    l_i = _l_i;
	    l_j = _l_j;
	    nnz = _nnz;
	}
	SpMatrix(int _l_i, int _l_j, int _nnz, Vector<T> &v) {
	    l_i = _l_i;
	    l_j = _l_j;
	    nnz = _nnz;
	}
	SpMatrix(const SpMatrix& spm) {

	}
	
	~SpMatrix() {
	    delete [] IA;
	    delete [] JA;
	    delete [] SA;
	    cout << "Have deleted." << endl;
	}
    };
    
    //get different Matrix
    template <typename T>
	Matrix<T> eye(int n) {
	Matrix<T> E(n, n);
	for (int i=0;i<n;i++) {
	    for (int j=0;j<n;j++) {
		if (i==j) {
		    E[i][i] = 1;
		} else {
		    E[i][j] = 0;
		}
	    }
	}
	return E;
    }
    
    template <typename T>
	Matrix<T> ones(int n) {
	Matrix<T> E(n, n);
	for (int i=0;i<n;i++) {
	    for (int j=0;j<n;j++) {
		E[i][j] = 1;
	    }
	}
	return E;
    }

    template <typename T>
	Matrix<T> zeros(int n) {
	Matrix<T> E(n, n);
	for (int i=0;i<n;i++) {
	    for (int j=0;j<n;j++) {
		E[i][j] = 0;
	    }
	}
	return E;
    }

    template <typename T>
	Matrix<T> hilbert(int n) {
	Matrix<T> A(n, n);
	for (int i=0;i<n;i++) {
	    for (int j=0;j<n;j++) {
		A[i][j] = 1.0 / (i+j+1);
	    }
	}
	return A;
    }

    template <typename T>
	Matrix<T> rand_mat(int row, int column) {
	Matrix<T> m(row, column);
	srand((int)time(0));	
	for (int i=0;i<row;i++) {
	    for (int j=0;j<column;j++) {
		m[i][j] = rand()/double(RAND_MAX);
	    }
	}
	return m;
    }
    
    template <typename T>    
    Vector<T> householder(Vector<T> &x) {
	int n = x.get_len();
        Vector<T> v(n);
	double model = x.inf_norm();
	
	//protect x
	if (model==0) {
	    return v;
	}

	x = x * (1.0/model);
	double sumOfx = x * x;
	double sumOfx1 = sumOfx - x[0] * x[0];
	v = x;
	v[0] = 0;
	double b;
	if (sumOfx1==0) {
	    b = 0;
	} else if (x[0]<=0) {
	    v[0] = x[0] - sqrt(sumOfx);
	} else {
	    v[0] = -1*sumOfx1/(x[0]+sqrt(sumOfx));	    
	}
	v = v / v.two_norm();
	return v;
    }    
}

#endif
